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ABSTRACT 


Several explanations of the vortex breakdown phenomenon have 
been proposed, including the finite transition theory of Benjamin 
and the boundary-layer separation analogy of Gartshore, Hall, and 
others. However, due essentially to a lack of quantitative in- 
formation with which to compare the predictions of these theories, 
none has gained general acceptance. A numerical investigation 
has, therefore, been undertaken to provide such information. 

Steady solutions of the Navier-Stokes equations, in terms of 
velocity and pressure, for breakdown in an unconfined viscous 
vortex are obtained numerically using the "artificial compressi- 
bility" technique of Chorin combined with an ADI finite-difference 
scheme, Axisymmetry is assumed and boundary conditions are care- 
fully applied at the boundaries of a large finite region in an 
axial plane while resolution near the axis is maintained by a 
coordinate transformation. The solutions, which are obtained 
for Reynolds numbers up to 200 based on the free-stream axial 
velocity and a characteristic core radius, show that breakdown 
results from the diffusion and convection of vorticity away from 
the vortex core which, because of the strong coupling between the 
circumferential and axial velocity fields in strongly swirling 
flows, can lead to stagnation and reversal of the axial flow near 
the axis. Breakdown is shown to be a necessary occurrence in 
flows with large enough values of a single parameter, calculated 
from conditions specified upstream, which decreases as Reynolds 
number increases toward the minimum required for the appearance 



of discontinuous solutions of the quasl-cylindrical equations, 
which apply only to high Reynolds number flows. These discon- 
tinuities are due to the large axial gradients which occur in 
flows which break down or nearly break down. Breakdowns in 
flows with both super- and subcritical upstream conditions are 
obtained so that the finite transition theory, which requires 
supercritical conditions upstream, is unable to fully account 
for the phenomenon. Based on the appearance in the solutions 
of a second retardation of the axial velocity behind the break- 
down bubble in highly swirling flows, an explanation of the 
several experimentally observed forms of breakdown is proposed. 
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1.0 INTRODUCTION 

1 .1 Motivation for This Study 

The rolled-up wake system behind a lifting wing consists 
essentially of two oppositely rotating viscous cores each of which 
is surrounded by a swirling Irrotational flow. In the case of the 
very heavy transport aircraft which have been introduced in the 
past several years, these viscous vortex systems can extend many 
miles, and remain vigorous, with high rotational velocities, for 
relatively long periods of time. They are a particular danger in 
the crowded airspace near airports, and encounters of smaller air- 
craft with such flows have resulted in numerous fatal accidents. 
There is an obvious necessity to in some way bring about their more 
rapid dissipation, and vortex breakdown is one of several mechanisms 
which have been suggested for this purpose. Unfortunately, the 
breakdown phenomenon is not very well understood, and while there 
is no lack of experimental observations, little quantitative 
information exists with which to correlate any of the several 
existing theories. 

It is with these considerations that a numerical approach 
to the study of vortex breakdown is undertaken. It is hoped that 
this study will yield practical information regarding the behavior 
of vortex flows exhibiting breakdown, and provide some insight 
into a poorly understood phenomenon. 

1 -2 Definition and Description of Vortex Breakdown 

To begin this study, a definition of vortex breakdown is 


required. 
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Vortex breakdown Is a pronounced decrease in the axial 
velocity on the axis of an axisymmetric swirling flow, and the 
corresponding divergence of the stream surfaces near the axis. 

Generally, the axial flow retardation results in stagna- 
tion and a region of reversed flow which may or may not be axi- 
symmetric. 

A summary of the experimental observations of Harvey 
(1962) and Sarpkaya (1971a, b) will serve to present the several 
forms in which vortex breakdown can appear. 

For flow in a tube, initially without swirl, a sequence 
of phenomena, made visible by injecting dye or smoke into the tube 
axis, may be observed as the amount of initial swirl is increased. 

The first to appear, for sufficient swirl, is an asymmet- 
ric phenomenon. The axial filament deforms into a spiral following 
an abrupt kink which occurs just upstream of a stagnation point. 

The spiral persists for a few turns, then breaks up into turbulence. 

As the swirl is increased a second asymmetric phenomenon 
may occur. The axial filament decelerates and expands into a 
curved sheet which deforms sinuously, and, as in the previous case, 
gradually breaks up into turbulence. There is no evidence of 
stagnation in these flows. 

For larger amounts of swirl, the axial filament spreads 
symmetrically at a stagnation point and the outer stream surfaces 
expand as if meeting a solid body. Behind the stagnation point 
there can be either a bubble of recirculating fluid followed by a 
return to conditions similar to those upstream and a spiral break- 
down like that described above, or a bubble immediately followed 
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by an expanding region of turbulent flow. 

Finally, if enough swirl is imparted, a core of reversed 
flow stretching along the axis for nearly the entire tube length 
can be obtained. 

1 .3 Observations of Vortex Breakdown 

Besides its possible occurrence in trailing vortex sys- 
tems, vortex breakdown has been observed and studied in a variety 
of flow situations. 

The first recorded observation of breakdown was in coni- 
cal leading edge vortices in separated flow over highly swept 
leading edges at high angles of attack (Peckham and Atkinson, 1957), 
and most of the early experimental investigations dealt with its 
effect on lift and moment coefficients (Elle, 1958; Werl^, 1960; 
Lambourne and Bryer, 1961), It was observed in supersonic as well 
as subsonic flow by Lambourne and Bryer (1961). 

Breakdown can occur In swirling flow in nozzles and 
diffusers (Gore and Ranz, 1964; So, 1967), and may also occur in 
tornado funnels (Morton, 1966), It was studied extensively by 
both Harvey (1962) and Sarpkaya (1971a, b) in vortex tubes with the 
aid of smoke and dye, as mentioned in the previous section. 

Although there is no lack of experimental observations 
of breakdown, the extreme sensitivity of the phenomenon to dis- 
turbances makes the introduction of probes impossible so that very 
little quantitative information exists. Most of what is known 
about flow conditions immediately upstream and downstream of break- 
down has been deduced from photographs. Laser anemometry, which 
is now being introduced to study trailing vortices, would be a 
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useful technique for studying breakdown. 

1.4 Theoretical Explanations 

The several theoretical explanations of vortex breakdown, 
which are discussed in detail in the review by Hall (1972), may be 
divided into three basic categories: 

(1) Vortex breakdown is similar to the separation of a 
two-dimensional boundary layer. 

(2) Vortex breakdown is a consequence of hydrodynamic 
instability, 

(3) Vortex breakdown is in some way dependent on the 
existence of a critical state. 

The first of these explanations, that breakdown is 
analogous to the separation of a two-dimensional boundary layer, 
developed from the attempts of Gartshore (1962,1963) and Hall 
(1966,1967) to obtain solutions of a system of equations which 
were obtained from the full, viscous, incompressible equations 
from the assumptions that the Reynolds number of the flow, based 
on a free stream velocity and vortex core radius, is large and 
that gradients in the axial direction are much smaller than in the 
radial. The resulting equations of this so-called "quasi-cylindri- 
cal" approximation are very much like the two-dimensional boundary- 
layer equations, and may be solved from given upstream conditions. 
Both Gartshore, who used an integral method similar to the Karmdn- 
Polhausen technique for the boundary layer equations, and Hall, 
who used a finite-difference method, found that for certain up- 
stream conditions their solutions could not be continued beyond 
some point downstream. They then assumed that this inability to 
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continue was due to the violation of the requirement that 
3/3z « 3/8r, and that breakdown must occur nearby. This is analo- 
gous to the assumption that boundary-layer separation occurs near 
the computed point at which wall shear vanishes even though the 
boundary-layer approximation may have failed somewhere upstream. 

The second theory, that breakdown is the result of a 
hydrodynamic instability to spiral disturbances was proposed by 
Ludwieg (1962,1965). He found a stability boundary for spiralling 
flows in a narrow annulus which was experimentally verified and 
suggested that it could be at least an approximate necessary con- 
dition for the stability of vortex flows. 

The third theory, that breakdown depends in some way on 
the existence of a critical state, was initially proposed by 
Squire (1960). He looked for conditions, which were later defined 
by Benjamin (1962) as the critical state, in which a standing 
cylindrical wave of infinite length could be first sustained, and 
assumed that such conditions could lead to breakdown. He con- 
sidered three different swirl distributions and found critical 
values of the ratio of the swirl velocity to the axial velocity 
to be of the order of unity. 

Benjamin (1962, 1965, 1967;Fraenkel ,1967), proposed that 
breakdown was a transition between two steady states of axisym- 
metric swirling flow and made an analogy with the hydraulic jump 
in open-channel flow. He showed that in an inviscid cylindrical 
swirling flow, the equation of motion for the stream function 
4) possesses an infinite number of solutions for given values 
of i|) at r = 0 and r = R, and given distributions. 
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with respect to i/j, of the circulation and the total head. He 
argued that only two of the solutions, which do not intersect 
except at their endpoints, are relevant, and showed that one of 
the solutions is supercritical, that is. cannot support infinites- 
imal standing waves, and that the other is subcritical, and can. 
The momentum flux (or flow force) 

R P 

2ir / (p + p w^) r dr 
0 

of the subcritical flow, however, is larger than that of the super- 
critical so that if a transition from one to the other is to occur 
the differences in the momentum flux must be made up. He then 
explained vortex breakdown as the transition of the flow from the 
supercritical to the subcritical state, with the appearance of 
axisymmetric standing waves in the subcritical state to make up 
the momentum flux difference when it is small. Larger momentum 
differences, however, cannot be made up in this way and he proposed 
that the leading wave then breaks and results in turbulence with a 
corresponding reduction in the momentum flux. 

Recent theoretical work has combined elements of these 

theories. 

Hall (1972) by considering two solutions, one expressed 
as a perturbation of the other, of the quasi-cylindrical equation 
for inviscid flow at cross sections a small distance iz apart, 
showed that the condition that the perturbation not vanish as Az 
approaches zero, that is, that axial gradients become infinitely 
large, is identical with the condition of the critical state. 
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He then described breakdown as the result of the initially super- 
critical flow, which may be described by the quasi-cylindrical 
approximation, approaching the critical state where the approxi- 
mation fails and derivatives increase without bounds. He thus 
identified with vortex breakdown both the failure of the quasi- 
cylindrical approximation and the appearance of the critical state. 

Mager (1972), following the work of Gartshore (1963) and 
Gossel (1967,1971), showed that the solutions of the quasi-cylindri- 
cal momentum-integral equations have two distinct branches and, 
with the appropriate upstream conditions, discontinuities which 
represent the appearance of infinite axial gradients and are 
assumed to signal the inception of vortex breakdown. He then 
demonstrated that Benjamin's finite transition could be explained 
as a sudden cross-over, upstream of the discontinuity, from one 
branch of the solution to the other, and suggested that the physi- 
cal manifestation of the discontinuity might be the spiral break- 
down while the cross-over, if it occurred, might appear as the 
axisymmetric bubble upstream. An Important aspect of this analysis 
is that the discontinuities may arise both in initially super- 
critical and subcritical flows, and that they occur in the vicinity 
of the critical state. 

The work of Hall and Mager suggests, therefore, that the 
critical state, which is determined from an inviscid analysis, is 
approximately the condition for which the quasi-cylindrical approxi- 
mation fails and breakdown occurs, and that this condition may be 
approached by both supercritical and subcritical flows. 
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1 .5 Outline of the Study 

This study obtains numerical solutions of the full 
Navier-Stokes equations for a single, unconfined, viscous vortex, 
imbedded in an irrotational flow with a uniform free-stream axial 
velocity. Chapter 2 describes the mathematical formulation of 
the problem, and Chapter 3 presents the numerical model. Results 
and discussion are presented in Chapter 4, and Chapter 5 summar- 
izes the findings of this study. 
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2.0 MATHEMATICAL MODEL 
2 J Equations 

While energy and angular momentum considerations require 
that an unconfined vortex appear as a member of an oppositely 
rotating vortex pair (Morton, 1969), a single vortex may be studied 
with the assumption that its opposite is far enough away to have 
little effect. Assuming axisymmetry- and that the flow is steady 
in a coordinate system fixed with respect to the vortex generator 
(e.g,, a lifting wing), the momentum equations for an incompressible 
flow In the radial, circumferential, and axial directions reduce to 


2 2 2 
rs(n M + w V 1 - lE. + S u 4. 9u . a u u > 


3z 


ar 


nfu + u iv . uy_ V . ^ / a^v . 1 av . a^v v . 


az 


p(u|^+wi^) = 1^.4) 




57 ■ " 3z ' 3z 
while mass conservation may be expressed as 


( 2 - 1 ) 

( 2 - 2 ) 

(2-3) 


i 3K) + iw , Q 
r ar az ^ 


(2-4) 


p is the density, y the viscosity, p the pressure, and u, v. 
and w are the velocity components in the r, 6, and z direc- 
tions, respectively (Fig. 1), 

Non-dimensional izing lengths by a characteristic core radius, 

★ 

r^,, velocities by W^, the free stream axial velocity, and defin- 
ing a dimensionless pressure, (p-p„)/pW^^, where is the 
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uniform ambient pressure infinitely far from the vortex, casts 
equations (2-1) through (2-4) in the following dimensionless form 


n^+w — - — = - + i ^ i* \ f?-5l 


Sv . 3v . uv _ 1 / 9^v .1 &v . V V 


dz 


„ 3w . aw _ ^ . 1 , a^ . 1 3w . a^w V 


(2-6) 


( 2 - 7 ) 


i 1Lmi1+ M = 0 

r 3r 3z “ 


( 2 - 8 ) 


where for convenience the same symbols for the variables have been 

★ 

retained, and Re = W^r^/v is the core Reynolds number. 

Although a solution of these equations is to be obtained for 
an unconfined vortex, boundary conditions applied at the boundary 
of a finite domain will be used for reasons discussed in a later 
section. 


2.2 Boundary Conditions 

Choosing the core radius at z = 0 as r* conditions 
are specified on the boundary of a region 0<r£R, 0<^z<L 
where R and L are much larger than one. The boundary condi- 
tions, depicted graphically in Fig. 2, are given below. A 
detailed explanation follows. 

At the upstream, z = 0, boundary the velocities are 
specified functions of r with 


for 


0 < r < R 


and 


u(r) = 0 


(2-9) 
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v(r) = V r(2-r^) 

for 0 <_ r £ 1 

w(r) = a + (1-a) (6 - 8r + 3r^) 

while 

v(r) = V/R 

for 1 < r < R 

w(r) = 1 

and a and V are defined below. 

At the downstream, z = L» boundary 

M - n 3v _ n _ n 

3^-0 , 3z ‘ ^ ^ ° 

are imposed for 0 < r £ R. 

At the axis, r = 0, 

u = 0 , v = 0 , and— =0 

for 0 £ z £ L 

while at the radial boundary, r = R 

=0 . V = V/R , and w = 1 (2-14) 

for 0 £ 2 £ L. 

The upstream, z = 0, conditions were chosen to approxi- 
mate the experimentally observed behavior of the velocity in the 
cores of trailing vortices (McCormick, Tangier, and Sherrieb, 1968; 
Chigier and Corsiglla, 1971,1972), and were used by Mager (1971) in 
his breakdown analysis. The vorticity is concentrated within the 
vortex core and the external angular velocity is then required to 
have the irrotational, r"\ distribution. V is the specified 


( 2 - 10 ) 


( 2 - 11 ) 


( 2 - 12 ) 


(2-13) 
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velocity at the core edge and is equal to the circulation around 
the core after non-dimensional ization by Zirr^W . The cubic form 
of v(r) allows solid-body-like rotation near the core center, 
and, since the derivatives are equal at r = 1, a smooth transi- 
tion to the irrotational flow at the core edge. The maximum 
angular velocity, 1,088 V, occurs at r = /273 , 

The quartic distribution of axial velocity in the core satis- 
fies the requirements that It join smoothly with the uniform axial 
flow outside, and that 9w/9r and, hence, shear, vanish at the 
axis. The form factor a is the ratio of the axial velocity at 
the core center to the velocity in the free stream. 

These distributions, requiring the specification of only two 
parameters, are quite flexible. Any desired amount of swirl and, 
therefore, circulation may be obtained by proper choice of V, 
while setting a less than or greater than one yields axial 
velocity profiles with axial momentum flux deficits or excesses 
(i.e., wake-like or jet-like profiles). Uniform axial flow 
results when a = 1 . 

The downstream, z = L, conditions assume that L has been 
chosen large enough that the rate of change of the flow properties 
with z is negligibly small. This approach is far less restric- 
tive than, for example, attempting to specify a priori the solution 
downstream. Although some small error is introduced, its effect 
upstream is limited for large enough L, as a numerical investi- 
gation described in a later section demonstrated. Similar con- 
ditions for primitive (v,p) variable formulations at outflow 
boundaries have been used by several authors (e.g., Taylor and 
Ndefo, 1970; Fortin, Peyret, Temam, 1971). 
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The axis, r = 0, conditions result from mass conserva- 
tion, which requires that u = 0, and from the requirement that 
shear stress at the axis and, therefore, v and 9w/3r vanish. 

Conditions at r = R are based on the assumption that 
phenomena near the axis, including both breakdown and diffusion of 
the vortex core have negligible effect on the circulation and 
axial velocity at large radial distances from the axis. The axial 
velocity is then simply that in the uniform free stream, while the 
swirl velocity is given by the irrotational distribution. Con- 
tinuity yields 3(ur)/3r = 0, allowing entrainment of fluid from 
outside the solution domain. 



14 


3.0 NUMERICAL METHOD 

3.1 Artificial Compressibility 

After considering a number of numerical methods for 
obtaining solutions of equations (2-5) through (2-8), Chorin's 
"artificial compressibility" technique (Chorin, 1967) was 
selected. This method solves the primitive variable equations 
directly, avoiding the use of vorticity-stream function equations 
which, for axisymmetric flows with swirl, are quite complicated, 
and which, in general, tend to yield less accurate numerical solu- 
tions than the corresponding primitive variable equations (Orzag, 
1971,1974). Chorin's method is simple, efficient, and has a sig- 
nificant advantage over other primitive variable techniques in that 
the solution of a Poisson equation for pressure and, consequently, 
boundary conditions for the pressure, are not required. A brief 
explanation of the method follows. 

Beginning with the equations of motion for a steady. 
Incompressible flow, for which a solution is desired 


v-Vv = - 7p + |^V^v (3-1) 

V • V = 0 (3-2) 

an auxiliary system of equations 

||-+v-Vv = - Vp + ^V^v (3-3) 

5 + V • V = 0 (3-4) 


is introduced. This new system features a time-like dependence on 
a new variable t, and the replacement of the kinematic constraint 
(3-2) by (3-4), which permits a simple explicit variation of the 
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pressure field. These auxiliary equations may be solved numer- 
ically from essentially arbitrary initial conditions to a steady 
(i.e., t-independent) limit using any of a large variety of 
finite difference schemes developed for initial value problems. 
Then, since steady solutions of (3-3) and (3-4) necessarily satisfy 
(3-1) and (3-2), the resulting numerical solution satisfies the 
original differential equations to the spatial order of accuracy 
of the particular difference scheme used. The desired solution 
is thus obtained. The artificial compressibility, 6, plays a 
role similar to a relaxation parameter and vanishes from the steady 
solution. Its numerical value is chosen for most rapid convergence. 

Convergence of the auxiliary system of equations has been 
proven for Stokes' flow with v specified at boundaries (Chorin, 
op. cit.; Fortin, Temam, and Peyret, 1971), and the method has been 
applied with various difference schemes to a number of problems 
(Chorin, op. cit.; Fortin, et al.; Plows, 1968). 

Applying the method directly to this study, solutions of 
(2-5) through (2-8) should be obtained by the solution of the cor- 
responding auxiliary equations, which are simply (2-5) through 
(2-8) with the additional terms 9u/3t, 9v/9t, 3w/3t, and 
63p/3t, respectively, carried out to a steady state. Rather than 
attempting a finite difference solution of the auxiliary equations 
directly, however, a coordinate transformation, described in the 
next section, was first applied. 

3.2 Coordinate Transformation 

Since vortex breakdown is essentially a phenomenon of the 
vortex core, good resolution and a finely-spaced grid are required 
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in and near the core. Unfortunately, uniform grid spacing in the 
rather large region, 0 < r < R, 0 < z < L, is then not possible 
because of the overwhelming number of grid points which would be 
required. However, explicit In the boundary conditions 


ly. 

3z 


= 0 



^ - 
3z 


0 


at z = L and both explicit and implicit in the condition 
=0 , V = V/R . w = 1 

at r = R are the assumptions that derivatives in the axial direc- 
tion far enough downstream and in the radial direction outside the 
core become increasingly small with r and z, respectively. 
Therefore, an unequally spaced array of grid points with increas- 
ing mesh width as both r and z increase is sufficient, and, 

since it makes more efficient use of each grid point, definitely 

advantageous. 

The effect of a non-uniform mesh has been obtained in 
this study through the application of the coordinate transformation 

y = i Jin (1 + r/b) (3-5) 

X = 1 Jin (1 + z/d) (3-6) 

which, with proper choice of a, b, c, and d, effectively stretches 
the region near the axis radially and the region near the upstream 
boundary axially, while contracting the more distant regions. 

These transformations have been used in a number of problems 
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(Skoglund and Gay, 1968; Rimon and Cheng, 1969; Pao and Daugherty, 
1969), and have been applied In this study to map the region 

0<z<L onto 0 < y £ 1/2, 0 £ x <_ 1 . Numerical 
values of a, b, c, and d are determined in practice from R, L, 
the total number of grid points in each direction, and the number 
of points desired in any Increment of r and z originating at 
zero. 

Details regarding the application of these transforma- 
tions to the auxiliary equation have been left to Appendix A, The 
transformed auxiliary equations are 
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where f = f(y), g = g(y), s = s(x), t = t(x) are coefficients 
defined in Appendix A, and e = l/r(y). 
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The transformed boundary conditions are as follows. 

For X - 0, 0 £ y _< 1/2: 

u, V, and w are specified functions of r(y) 

as given previously. (3-^11) 

For X = 1 , 0 £ y £ 1/2: 

= 0 . = 0 , and = 0 . (3-12) 

For 0 £ X £ 1 , y = 0: 

u = 0 , v = 0 . and|j=0 (3-13) 

while, for 0 £ x £ 1 , y = 1/2: 

“ e ly ' ° • V = V/R and w = 1 (3-14) 

The finite difference scheme developed to solve equations (3-7) 
through (3-10) subject to the boundary conditions above is dis- 
cussed in a later section. 

3.3 Exploratory Calculations 

To obtain a numerical solution for vortex breakdown 
several approaches, briefly described in the following paragraphs, 
were explored. Each provided essential information and experience 
without which the final scheme used to obtain the solution presented 
later in this report could not have been developed. 

Initially an explicit finite difference scheme with 
simple programming was constructed. However, instead of (3-5) 
and (3-6), the transformation 



was applied. This was used to map the infinite region 0 < r < «>, 
0 z < » onto 0£y£l, 0£X£l so that conditions u = 0, 
V = 0, and w = 1, which apply as r and z approach infinity 
could be used at the boundaries of the transformed region without 
approximation. The form of the z to x transformation was 
chosen to match, in the sense that singularities are avoided, the 
asymptotic solution of Batchelor (1964) for a viscous trailing 
vortex. The r to y transformation was chosen for its easy 
inversion. The resulting transformed auxiliary equation, which, 
except for the definitions of the coefficients f, g, etc., are 
almost identical to (3-7) through (3-10), were solved by an 
explicit scheme which used centered differences for the diffusion 
and pressure terms, forward differencing for t, and first-order 
upwind differencing for the convection terms. These Initial coarse 
grid calculations established the possibility of obtaining and 
studying a vortex breakdown numerically at Reynolds numbers 
sufficiently low for finite difference solutions to be meaningful. 

In order that solutions might be obtained more effici- 
ently other finite difference methods were considered and an 
implicit scheme based on the alternating-direction implicit (ADI) 
method of Peaceman and Rachford (1955) was eventually selected. 

This scheme was significantly faster. Time steps much 
larger than possible with the explicit scheme could be used, easily 
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compensating for the extra computational labor involved. Initially 
all spatial derivatives were approximated by centered differences, 
but experiments, even with a fine mesh spacing, 61 points in the 
x-direction and 31 in the y-direction, demonstrated that only 
solutions for very low Reynolds numbers (< 50) could be obtained. 
However, by applying upwind differencing for x > 2/3, flows at 
Re = 100 could be calculated. Higher Reynolds number flows re- 
quired upwind differencing everywhere. Fully aware that there 
could be significant errors due to the large numerical viscosity 
introduced by upwind differencing, calculations were performed for 
Reynolds numbers from 100 to 1000. These computations later proved 
essential in providing estimates for the approximate location of 
breakdown, the breakdown bubble's radius and length, and the magni- 
tude of the reversed axial flow in the bubble. 

The reason only extremely low Reynolds number flows 
could be obtained without upwind differencing is apparent. While 
the grid spacing in the transformed (x,y) plane could be made 
arbitrarily small through the use of an increased number of grid 
points, the grid spacing as it would appear in the r-z plane 
always becomes very large and eventually infinite as x and y 
approach one. Letting represent the mesh spacing in the r-z 
plane, the criterion 

Re = 0(1) 

which has been claimed to be a necessary condition for a numerical 
solution (Dennis and Cheng, 1969; Chorin, 1972), will always be 
violated for large r and z. A similar criterion for a solution 
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In x-y space Is obtained by letting h represent the fixed grid 
spacing in that space. Then 

h _ Ax ^ dx 

k;: - u n 

and since dx/dz = s(x), the coefficient of the transformed 
equation defined in Chapter 2, 



Therefore, the condition 

Rei = 0(1) 

must be satisfied. As x approaches one, s approaches zero 
and this criterion is violated. For the Re - 50 case mentioned 
above, the violation occurs at a small number of downstream grid 
locations, and apparently a solution could still be obtained. 

For Re = 100, however, violation was more widespread, and a solu- 
tion could only be obtained by using upwind differencing for 
x > 2/3 which, in effect, lowered the Reynolds number in that 
region with the introduction of numerical viscosity. 

An examination of the solutions obtained in these calcu- 
lations suggested that approximations could be made at finite r 
and 2 . Boundary conditions could then be formulated for a finite 
region in r-z space, while allowing a solution to be obtained 
for an unconfined flow. Use of such a finite domain was expected 
to avoid the mesh width difficulties explained above, and to allow 
the calculation of accurate solutions at higher Reynolds numbers 
than possible in the Infinite domain. The finite domain model of 
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Section 2.2 and the finite transformation presented in 3.2 were 
then developed. 

3.4 The Finite Difference Scheme 

An ADI scheme similar to that described in the previous 
section for the infinite domain calculations was used to obtain 
steady solutions of (3=7), (3-8), (3-9) and (3-10), subject to 
the boundary conditions (3-11) through (3-14), for a wide range 
of the parameters Re, a, and V. The scheme, including treat- 
ment of the boundary conditions, is presented in detail in Appen- 
dix B. Centered differences were used for all spatial derivatives 
Including convection terms, so that no numerical viscosity was 
introduced. A staggered grid system, which was used to increase 
the accuracy of the scheme and to avoid the necessity of calculat- 
ing pressure at the boundaries, is also described in the appendix. 
The spatial truncation error of the scheme is O(Ax^) + O(Ay^) and, 
except for the effects of first order approximations for the 
derivative boundary conditions, solutions are of the same order. 

The location of the radial and downstream boundaries, 
r = R and 2 = L, were determined in the following manner. 
Estimates were first obtained from an examination of the results 
of the infinite domain calcaulations at Re = 100. Then, using a 
grid system with 61 points in the x-direction, 31 points in the 
y-direction, and 12 points in the intervals 0<r£l,0<^z<l, 
solutions with various values of R and L, based on the esti- 
mated values, were obtained. These solutions demonstrated that 
R could be safely set at 10 and L at 20 since the use of larger 
values did not significantly change the solutions, and, also. 
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since a comparison of these solutions with those obtained in the 
infinite domain showed very good agreement. Such a comparison is 
presented in Fig, 4 for the axial velocity on the axis versus z. 

Numerical experiments which varied 6 in the hope of 
obtaining an optimum value demonstrated that larger values of 6 
permitted longer time steps, while shorter steps were required 
for smaller This behavior was predicted by Chorin. Very long 
time steps based on large values of 6 did not permit stable 
calculations to convergence, probably due to the timewise lineari- 
zation of the convection terms in the finite difference momentum 
equations. Eventually, 6 was set equal to one, and a time step 
of length 0.11 was then used in all calculations. 

Computations were run to convergence, defined as changes 
in velocities less than 0,0005 over 100 time steps with the 61 x 31 
grid system explained above. This required from two to four 
thousand time steps, depending on the particular case and initial 
conditions. This rather large number of time steps resulted from 
the quite complicated three dimensional dynamics of the flow in 
vortex breakdown. When possible, the initial condition of any 
computation was a solution from a previous computation with 
values of Re, a, and V close to those of the case being con- 
sidered. This was most efficient, saving as many as 1500 time 
steps, for fixed Re and a, as V was varied from low to high 
values. Each computation required from two to four minutes of 
computer time on the Lawrence Berkeley Laboratory CDC 7600 when 
compiled with the CDC FTN optimizing compiler. Graphical output 
was disposed to both microfiche and 35mm film, from which many of 
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the figures presented later in this report were obtained. 
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4.0 RESULTS AND DISCUSSION 

4.1 Preliminary Comments 

A total of 39 solutions with various combinations of Re, 
a, and V were obtained. The combinations, which are summarized 
in Table 1, were chosen with the Intent that each solution should 
be as informative as possible, and that the net effect of all the 
solutions should be to provide a comparison of the several existing 
theories of vortex breakdown with quantitative results. 

It should be noted prior to a discussion of the solutions 
that centered differencing of all spatial derivatives was used for 
computations with Reynolds numbers up to 200. Solutions at higher 
Reynolds numbers, which were obtained for qualitative purposes 
only, required upwind differencing for the convection terms in the 
momentum equations and, therefore, a slight modification of the 
program described in the previous chapter. 

4.2 Solutions for Vortex Flow 

The development of a vortex breakdown, in terms of a 
sequence of steady state solutions with increasing amounts of swirl 
will be described in this section. These solutions were all 
obtained with Re = 200, a = 1 so that the axial velocity dis- 
tribution at z = 0 was uniform. Calculations were performed 
with V set at 0.63, 0.80, 0,85, 0.894, 1.0, and 1.095. In most 
cases, little of interest occurs outside the region 0 f_r ^2, 

0 £ z £ 6, which is a small part of the domain 0 £ r £ 10, 

0 120, in which the solutions have been obtained. Therefore, 

most of the figures to be presented will display, in detail, various 
aspects of the solutions in this small region only. 
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Figure 5a shows a contour plot for a number of constant 
values of the Stokes stream function, i|), defined by 

(r,z) 

= / r(w dz - u dr) 

where the line integral is taken along an arbitrary curve in an 
axial plane joining a reference point where = 0 to the point 
(r,z) (Batchelor, 1970). was computed at each grid point by 
simple numerical integration of the velocities from the solution 
obtained with V = 0.63. Each contour represents a stream surface 
and, since the flow is steady, a fluid particle entering the domain 
at z = 0 proceeds downstream spiralling along the stream surface 
on which it enters. Plots of other aspects of this solution are 
presented in Figs, 5b through 5g, 

Figure 5b displays axial velocity profiles at a number 
of axial locations, and Fig. 5c shows the variation of the axial 
velocity at several values of r with distance downstream. The 
slight retardation of the axial velocity within the core reduces 
the axial velocity at r = 0 to about 0.80 at z = 6. 

Figure 5d shows swirl velocity profiles at the same 
axial locations chosen for Fig. 5b, and Fig. 5e displays the 
variation of the swirl velocity with z at the same radial posi- 
tions plotted in Fig. 5c. The swirl velocity at r = 0 is always 
zero and not labeled in Fig. 5e. A reduction in the swirl veloci- 
ties in the core due to the diffusion of vorticity away from the 
axis is evident in these figures. For example, the velocity at 
r = 0.486 decreases from 0,55 at z = 0 to 0.48 at z = 6. 

It should be noted that the changing velocity field, as 
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Morton (1969) has pointed out, associated with vorticity diffusion, 
is due to a stress distribution which, in effect, transmits inward 
a retarding torque resulting from the non-zero circulation at radial 
infinity. The angular momentum of a disk of fluid of unit thick- 
ness, concentric with the axis, will, therefore, decrease with z. 
This may be clearly demonstrated by integrating the angular momen- 
turn equation (2-2), from r = 0 to r = r , where r > r , to 
obtain, after rearrangement 

* 

bI / It ^ *'5 * I" - 1? 

u r 

where r is the fixed constant circulation outside the core. The 
angular momentum flux along a cylinder of radius r* is, there- 
fore, not fixed but increases due to entrainment of fluid from the 

exterior of the cylinder, and decreases due to the effect of a 
* 

torque at r , and, although the angular momentum of the disk 
becomes infinite as r approaches Infinity, it must still de- 
crease at a finite rate with z since (ur) * does not, in 

y =00 

general, vanish. The flux of angular momentum, therefore, is not 
constant In vortex flows. 

Figure 5f presents the variation with z of three 

quantities: the pressure, p^, along the row of staggered grid 

points nearest the axis; the radius of the vortex core, r , 

c* 

which has been defined here as the cylindrical region, concentric 
with the axis around which the circulation is 0.95 V; and the 
maximum swirl velocity in the core, Note that since the 

pressure within the core will be less than the ambient pressure. 
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2 

the non-dimensional pressure (p-pj/pw^ will always be less 
than zero. In this case p^ varies very slowly with 2 , increas- 
ing from about -0.63 at 2 = 0 to -0.50 at z = 6, Diffusion 
of vorticity away from the axis is responsible for both the increase 
of the core radius from about 0.90 at z = 0 to 1.20 at z = 6, 
and the reduction of the maximum swirl velocity from 0,69 to about 
0.57. 

The axial velocity decrease within the core, noted above, 
may now be explained. As vorticity diffuses beyond the initial 
core radius, the swirl velocity of a particle of fluid near the 
axis will decrease. The radial pressure gradient required to 
balance the centrifugal acceleration of the particle must, there- 
fore, also decrease. This reduction is obtained through the 
increase in the pressure near the center of the core, displayed in 
Fig, 5f. The axial flow near the axis, encountering a small ad- 
verse gradient of pressure is, therefore, slightly retarded. 

The above argument may be demonstrated analytically 
using the approximation 



which expresses the balance of a particle's centrifugal accelera- 
tion with the restraining pressure force. This equation is exact 
only for slowly expanding cores in the limit of high Reynolds 
numbers, and forms part of the "quasi-cylindrical" approximation 
of the Navier-Stokes equations which is valid for such conditions 
(Gartshore, 1963; Hall, 1966; Bossel, 1967). Integration with 
respect to r from the axis to radial infinity followed by 
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differentiation with respect to z, yields 



-P/ f 


0 


where is the pressure at the center of the core. Since v 
is always positive, and viscous dissipation of swirl makes 3v/3z<0 
there results a positive axial pressure gradient at the axis, and 
thus the retardation of the axial flow near the axis. 

Figure 5g shows the axial velocity at the axis for the 
entire length of the computational domain. The deceleration of 
the axial flow does not continue indefinitely. As v and -3v/3z 
decrease, the rate of change of the pressure will also decrease so 
that the adverse pressure gradient becomes smaller with z. This 
behavior is predicted by equation (4-2). Viscous forces are then 
able to overcome the pressure gradient and accelerate the fluid, 
eventually to the free stream axial velocity. 

The solution obtained with V = 0.80 is displayed in 
Figs. 6a through 6g. 

A slight expansion of the stream surfaces in Fig. 6a 
is barely apparent so that the solution seems little different 
from that obtained with V = 0.63. However, the deceleration of 
the axial flow in the core (Figs, 6b and c) is much more pronounced 
with the velocity at r = 0, for example, decreasing to 0.55 at 
z = 6. Also more pronounced is the rate at which the swirl 
velocities (Figs. 6d and e) decrease with z. The swirl at 
r = 0.486 decreases from 0.70 at z = 0 to 0.41 at z = 6. 

This increased rate of swirl velocity decrease is due to the 
increased concentration of vorticity near the axis, which results 



30 


results In more rapid diffusion. In terms of the shear stress 

2 

distribution, t at z = 0 is equal to -2yVr so that the 

torque applied to a disk of fluid is of unit thickness, for r £l, 
4 

-47naVr , directly proportional to V. 

Using, once again, the notion of a balance between the 
centrifugal acceleration of fluid particles and the radial pressure 
gradient, the increased rate of swirl reduction with z must be 
supported by a more rapidly increasing pressure on the axis. Since 
V and -3v/9z are larger than in the previous case, 3 Pq/3z must 
also be larger. Comparison of 6f with 5f shows the difference in 
the pressure variation. For V = 0.63 the pressure varies only 
very slowly with z, while for V = 0.80 the pressure, which 
initially must be lower to provide the larger radial pressure 
gradient required for the increased swirl, shows a much more rapid 
rise, increasing from -1,03 at z = 0 to -0.63 at z = 6, The 
more pronounced retardation of the axial velocity noted above 
corresponds to this larger adverse pressure gradient. 

Figure 6f also displays the increasing core radius 
and the decreasing maximum swirl velocity in the core. Figure 6g 
shows that the axial velocity decreases to a minimum of about 0.5 
before shear forces begin the acceleration which will eventually 
result in a uniform axial flow Infinitely far downstream. 

Increasing V from 0,80 to 0.85 results In the solution 
presented in Figs. 7a through 7g. The stream surface expansion 
(Fig. 7a) is now clearly evident. The radii of the surfaces near 
the axis reach a maximum at about z = 5, and the axial velocity 
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on the axis (Figs. 7b and c) is quickly decelerated to a minimum of 
0.25 also near z = 5. The swirl velocity in the core (Figs. 7d 
and e) after decreasing initially in a manner very similar to the 
previous cases, reaches a minimum near z = 5, where, for example, 
the velocity at r = 0.486 is 0.28, and then begins to increase. 

The acceleration of the swirl velocity near the center of the core 
continues for a short distance, then, although not shown in Fig. 

7e, gives way to a final deceleration, which eventually will lead 
to total dissipation. Figure 7f shows that after a rapid initial 
increase, the pressure reaches a constant value, -0.57, and 
begins to decrease slightly, also near z = 5. 

The dynamic situation may be explained as follows. A 
fluid particle enters the domain at z = 0 close to the axis with 
axial velocity and angular momentum. As it proceeds downstream, 
shear forces reduce its angular momentum and, as explained previ- 
ously, this is associated with an adverse gradient of pressure which 
acts to decrease the axial velocity of the particle. However, since 
both V and -3v/3z are initially larger than the previous cases, 
the amount of retardation of the axial flow is also larger, and 
continuity requires that a significant amount of fluid move radi- 
ally outward, This movement is apparent in the stream surface plot. 
Then, since, except for viscous effects the angular momentum 
of the fluid particle is conserved as it moves radially, its angu- 
lar velocity must decrease. Figure 7d illustrates this behavior. 

As z increases from zero to about 5, the swirl velocity at any 
radial position is that of a fluid particle which has originated 
at some position closer to the axis. The effect of this is the 



32 


development of a slowly rotating region on the axis apparent from 
about 2=4, For 2 less than about 3, the reduction of swirl 
due to the radial outflow explained above requires, when combined 
with viscous effects, a much larger increase in the pressure at 
the axis. This is clearly seen in Figure 7f, For larger 2 , 
however, since v and 3v/32 are very small near the axis, the 
adverse pressure gradient also becomes small, and shear forces are 
then able to accelerate the axial flow near the center of the core 
(Fig, 7c). Continuity then requires an inward motion of fluid 
toward the axis, and, again assuming angular momentum conservation 
for a particle, the increase in the angular velocity at radial 
psoitions, displayed on Fig. 7e, must occur. The increasing 
angular velocity near the axis requires that the pressure on the 
axis decrease slightly, as may be seen in Fig. 7f at z = 6. 

Beyond this point the pressure continues to slowly decrease until 
approximately 2 = 7.5, assisting the shear forces in accelerating 
the axial flow (see Fig. 7g), Then the effect of viscous diffusion 
once again becomes significant, reducing the swirl velocities near 
the axis and leading to an adverse pressure gradient which decreases 
the rate of acceleration of the axial velocity. This is apparent 
in Fig. 7g from about 2 = 7 to 8. As in previous cases, however, 
this adverse gradient decreases with z as v and -3v/32 de- 
crease so that the axial velocity will eventually be accelerated 
by shear forces to the free stream velocity. 

The solution obtained with V = 0.8944 is presented in 
Figs. 8a through 8g. A small, well-developed bulge in the stream 
surface near 2 = 3 is readily apparent in Fig. 8a, Similar stream 
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surface expansions were obtained experimentally by Sarpkaya (1971), 
and analytically by Bossel (1967) for flow conditions close to 
those for which his integral analysis of the quasi-cylindrical 
equation failed. 

The axial velocity devrease (Figs. 8b and c) in the core 
is considerably more rapid than those enouuntered in the previous 
cases, and nearly results in stagnation at 2=3 where the axial 
velocity on the axis reaches a minimum of 0.02. The acceleration, 
then deceleration of fluid outside the core as it passes over the 
bulge may be detected in Fig. 8b, and is seen clearly at r = 1 
in Fig. 8c. 

The swirl velocity profiles (Figs. 8d and e) show, as 
in the previous case, the decrease in the swirl velocity due to 
shear stress and, especially, the reduction induced by the stream 
surface expansion. The reduction near the axis is particularly 
substantial. As Fig. 8c demonstrates, the swirl velocity at 
r = 0.486 decreases from 0.77 at z = 0 to 0.20 at z = 3. 

As shown in Fig. 8f, the pressure increases steeply from 
-1.21 at 2=0 to -0.76 at z = 2, and then becomes essentially 
constant within the bulge region, 2 < z £ 4, with a value of 
about -0.65. The contraction of the stream surfaces, begun by 
shear forces near the axis, as explained previously, then causes 
the pressure to decrease to -0.72 at z = 5.5 as the angular 
velocity of fluid near the axis increases. Viscous effects, how- 
ever, begin to decrease the swirl velocity and, once again, an 
increasing pressure near the axis is required. 

Also displayed in Fig. 8f is the rapid increase in r 
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between z = 2 and 3, which is due in the most part to radial con- 
vection of vorticity; andg the sharp decrease in from 0.96 

at 2 = 1 to 0.77 at z = 3. 

Figure 8g presents the variation of axial velocity on the 
axis. The Initial decrease, as explained above, is due to the large 
adverse pressure gradient, and the acceleration which begins at 
z = 2 is Initially due to shear forces but given assistance, 
beginning at about z = 4, by a favorable pressure gradient. The 
second deceleration, beginning at about z = 5.5, is the result 
of the second adverse pressure gradient and the following acceler- 
ation is due to the effect of shear forces as this gradient 
decreases. 

Figures 9a through 9g display a vortex breakdown obtained 
with V = 1.0, The small stream surface surface bulge of the pre- 
vious solution has developed into a closed bubble of recirculating 
fluid, similar in shape to some of those obtained by Sarpkaya (1971a), 
The stream surfaces within the bubble have small negative values. 

The outer stream surfaces, after expanding over the bubble, neck 
down behind it, and expand slightly again. 

The axial velocity plots (Figs. 9b and c) show the very 
rapid deceleration near the axis upstream of the bubble, w at 
r = 0 reaches a minimum of -0.082 at z = 1.64, then increases 
to 0.359 at 2 = 3,62, A second deceleration reduces the velocity 
to 0.176 at z = 5,16, where it begins its gradual return to free 
stream conditions (Fig. 9g), Clear in both Figs. 9b and 9c is the 
acceleration of fluid outside the core as it passes over the 
bubble. 
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Figures 9d and 9e denionstrate that due to the stream sur- 
face expansion the swirl velocities within and near the bubble are 
very small. At r = 0.486, for example, the swirl decreases from 
0.86 Initially to 0.124 at z = 1.78, before increasing to 0.53 at 
z = 3.96 due to the stream surface contraction behind the bubble. 

Figure 9f shows a very steep increase in the pressure 
near the axis ahead of the bubble. The pressure near the axis 
within the bubble is essentially constant, however, due to the very 
slow rotational speeds mentioned above. As in the previous cases, 
the stream surface contraction results in a pressure decrease, 
beginning at z = 1.78, which assists the acceleration of the 
axial velocity. Beyond z = 3.6, the pressure begins to increase 
with the dissipation of swirl due to shear stresses. 

The solution for a large, more vigorous breakdown, ob- 
tained with V * 1,095, Is displayed in Figs. 10a through lOg. 

The bubble, as shown in Fig. 10a, is larger than the 
previous case and has moved further upstream. The stream surface 
expansion around the bubble is now very rapid, and the contraction 
behind the bubble is followed by a second pronounced expansion, 
beginning at about z == 3. The near entrainment into the bubble 
of the stream surface closest to the axis should especially be 
noted. Sarpkaya (1971a) reported that fluid did not enter through 
the front of his experimentally obtained bubble, but instead 
passed over it and entered from the back. Then, after mixing 
turbulently inside, fluid exited into a second core-like region 
behind the bubble. While turbulent mixing has not, of course, been 
included in this study, the similarity between Fig. 8a and his 
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photographs is striking. 

Further comparison of Fig* Ba with Sarpkaya's photographs 
suggests that the second expansion might be the axisymmetric 
counterpart to the spiral breakdown which is always found, esiperi- 
mentally, behind an axisymmetrlc breakdown. If the flow in this 
second expansion were unstable to spiral disturbances, which of 
course do not exist in these numerical solutions, such a spiralling 
breakdown might then develop. This hypothesis will be discussed 
in somewhat more detail <in a later section. 

Figures 10b through lOg show that the solution differs 
from those previously discussed in some interesting respects. 

An increased rate of deceleration of the axial flow com- 
pared to the previous case is apparent in Figs. 10b and 10c 
although the minimum velocity attained, -0.036 is somewhat larger. 
Also, the axial velocity on the axis is reversed for a shorter 
distance than before. 

The swirl velocity near the axis (Figs. lOd and lOe) 
shows a more rapid decrease with z, which corresponds to the 
more rapid stream tube expansion. The bubble is clearly (Fig. lOd) 
a region of very slow rotation despite the very large initial 
swirl . 

Figure lOf displays the very large initial pressure 
gradient and, as in the previous case, the almost constant pressure 
in the bulge, which is followed by a decrease due to the stream 
surface contraction. Also, displayed in Fig. lOf are the core 
radius, which Increases rapidly for small z due to outward radial 
convection of vorticity, then decreases slightly due to inward 
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convection during the contraction, and which decreases 

iTldX 

quickly due to the stream surface expansion. 

Following the discussions of the next two sections, which 
deal with the solutions presented here, solutions for other 
values of a, V, and Re will be discussed, 

4.3 The Effect of Swirl 

It is clear from the solutions presented in the previous 
section that the appearance in a vortex core of a large adverse 
pressure gradient, with a corresponding large decrease in axial 
velocity, is critically dependent on the magnitude of the swirl 
which is introduced at z = 0. Figure 11 has been prepared to 
quantify this dependence. It shows the minimum axial velocity, 

W^, attained on the axis versus V for each of the solutions. 

The dashed trend curve represents a graphic extrapolation of the 
solutions to small values of V. As V increases from zero, where 
no axial retardation will occur, dW^^/dV becomes increasingly 
negative so that for small V, small changes in V yield only 
small changes in but for larger V, particularly near 

V = 0.80, small changes in V yield very large changes in 
as may be seen in the figure. An explanation of this highly non- 
linear behavior may be as follows. 

For small values of V, the effect of swirl dissipation 
is to increase the pressure on the axis and thus cause a small 
amount of axial flow retardation. The small radial flux away from 
the axis which is required by continuity, has, however, little 
little effect on the swirl velocity distribution. The swirl 
velocity field Is, therefore, essentially decoupled from the rest 
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of the flow field and affects it only through the small adverse 
pressure gradient. For larger amounts of swirl, however, this 
decoupling is no longer possible. The axial velocity retardation 
arising from viscous dissipation of the swirl is large enough to 
require a significant amount of radial outflow, and, as explained 
in the previous section, this decreases the swirl velocities near 
the axis, and thus supports a very much larger adverse pressure 
gradient and an increased reduction of the axial velocity. The 
very large negative value of dW^dV which occurs at about V = 0.8 
for the values of a and Re of these solutions is, therefore, due 
to the coupling, as V increases, of the initially decoupled swirl 
velocity field with the rest of the flow field. 

4.4 Comparisons with Experiments and Theory 

Before considering other solutions of the equations for 
breakdown, it is worth comparing the previously discussed solutions 
with experimental observation and the theoretical explanations 
described in Chapter 1. 

As pointed out earlier, there is very little experimental 
data regarding vortex breakdown. The quantitative information 
which does exist has, for the most part, been obtained from photo- 
graphs taken of flows, after the introduction of dye or smoke. In 
this way Harvey found that the swirl angle, <j> = tan“^(v/w), ahead 
of his breakdowns was never larger than 50.5®, and in a similar 
experiment Sarpkaya (1971a) found that it was never larger than 51®. 
Theoretical predictions of based on the finite transition con- 
cept, vary from 45® (Squire, 1960) to 62.5® (Bossel, 1967), depend- 
ing on the velocity profiles chosen for analysis. This is a 
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considerable range of variation, particularly since the maximum 
value of v/w, from which 4* is obtained, varies, therefore, 
from 1 .0 to 1 ,88, 

The maximum swirl angle at z = 0 of the solutions dis- 
cussed in Section 4,2 varied from 34.38° to 49.95° (see Table 2), 
and was 44,17° in the solution obtained with V = 0,8944 in which 
the axial flow almost stagnated. There is, therefore, reasonable 
agreement in terms of the swirl angle between the experimental 
evidence, theoretical predictions, and the numerical solutions. 

It should be noted that Mager (1972) substituted the 
Initial velocity profiles used in his analysis and in this study, 
into Benjamin's test equation for criticality and solved it with 
various combinations of a and V. He was then able to divide an 
a-V plane into supercritical and subcritical regimes. Based on 
his results, the upstream conditions for all the numerical solu- 
tions presented thus far are supercritical, as required by the 
finite transition theory for breakdown to occur. 

The second theoretical explanation reviewed in Chapter 1, 
that breakdown is the result of a hydrodynamic instability to a 
spiral disturbance, is obviously irrelevant since no such dis- 
turbances exist in these axisyiranetric calculations. As mentioned 
in Section 4,2 and also by Hall (1972), such an instability may 
determine, however, which type of breakdown, spiral or axisym- 
metric, occurs in experiment. 

The third theoretical approach to breakdown, that it is 
a swparation-like phenomenon, predicted by the failure of the quasi- 
cyllndrical approximation, must also be considered, Mager (1972) 
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provided a very precise criterion for the appearance of discontinu- 
ities, which he assumed were the results of such failures, in his 
integral solutions of the quasl-cylindrical equations. He defined 
a parameter 6-j which is the Invariant flux of axial momentum 
deficiency In the core plus a term dependent on the circulation 
around the core plus a third term whose value changes only behind 
the breakdown, Is determined from the upstream condition, and 
depends in no way on the super- or subcriticality of these condi- 
tions. In general increasing a decreases 6^ since the axial 
momentum flux is then larger, while increasing V increases 6-j, 
since the pressure deficit in the core is then larger. Mager then 
showed that if the Initial conditions, as determined from the values 

of a and V, are such that 6-j Is less than a certain value, 

★ 

, then there is a continuous solution. If, however, , 

discontinuities will arise, marked by the appearance of infinite 
gradients regardless of the sub- or supercriticality of the upstream 
flow. He determined that 6^* = -0.163989. 

Values of for all the solutions of Section 4,2, are 
shown in Table 2. From the discussion of Section 4.2 it is apparent 
that flow reversal in the numerical solutions corresponds to a value 
of 9^ slightly larger than -0.056, which is significantly larger 
than 0^ . There are a number of possible explanations for this 
difference. First, it is not clear that the failure of the quasi- 
cylindrical equations requires fully reversed flow, since such a 
failure represents only the inability of the approximation to deal 
with large axial gradients, and, as noted in Section 4„2, large 
gradients do appear in solutions which do not exhibit stagnation. 
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6^ could then possibly be a good deal less than the required 
for flow reversal, as is the case in these solutions. Second, 
Mager’s assumed profiles may not be general enough to describe the 
flow approaching breakdown. This is a possibility, as noted by 
Hall (1972), in all integral formulations. Third, and most obvious, 
Mager's result applies to high Reynolds number flows and 200 is, 
at best, only moderate. 

If 6-| is, however, the essential parameter, the occur- 
rence of breakdown should be independent, except for Reynolds num- 
ber effects, of the particular values chosen for a and V, and 
depend only on the value of . Also, if the theoretical explana- 
tion proposed in the first chapter, that breakdown is a result of 
initially supercritical or subcritical flows approaching a critical 
state in which very large axial gradients develop, is to be demon- 
strated, breakdowns in initially subcritical flows must be obtained. 
These ideas were important in the choice of which other solutions 
would be obtained, 

4,5 Other Solutions 

Several solutions for Re = 200 not described in Section 
4.2 are displayed in Figs. 12 through 18. These solutions are, in 
most respects, quite similar to those previously discussed. 

Figures 12 and 13 show the stream surfaces for two solu- 
tions with a = 0.3. Since a is very small, the axial momentum 
of the flow near the axis will also be small, so that the adverse 
pressure gradient which develops as a result of swirl dissipation 
need not be very large to produce stagnation. The solution with 
V = 0.8944, therefore, exhibits reversed flow (Fig. 13), while 
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that with the same swirl but at a = 1, discussed previously 
(Fig. 8) did not. It should also be noted that this solution 
results from upstream conditions which are subcritical. The 
appearance of breakdown in subcritical flows, as required by the 
proposed explanation discussed previously, is thus demonstrated. 

Solutions obtained with a = 0.6 and V = 0.82, 0,8944, 
and 1.095 are shown in Figures 14, 15, and 16. While the solution 
for V = 0.85 and a = 1 discussed previously showed only very 
slight stream surface expansion, the expansion is quite pronounced 
for V = 0.82 (Fig. 14), and as with a = 0.3, setting V to 
0.8944 is sufficient to produce reversed flow (Fig. 15). Figure 
16a displays stream surfaces obtained from the solution with 
V = 1.095 for which the upstream conditions were again subcriti- 
cal. The second retardation of the axial flow (Fig. 16b), which 
always appears in solutions with large swirl velocities, is, in 
this case, sufficient to produce stagnation and a second reversed 
flow region. 

Figures 17 and 18 show solutions obtained with a = 1.4, 
and it is apparent that V 1n excess of 0,95 is required for flow 
reversal. In this case, the axial momentum of the flow near the 
axis is very large and, therefore, a very strong adverse pressure 
gradient, which will exist only for large values of V, is 
required to produce stagnation. 

Attempts were made to obtain solutions for a greater 
than 1.4. These invariably failed for values of V large enough 
to be of Interest in this study. Also, solutions with V greater 
than about 1.2 for a = 1, which were desired to study further 
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flows with high swirl, failed. These failures were probably due 

to an incorrect scaling of the equations for flows with high swirl 

or large axial velocity excesses. For flows in which a or V 

are large, the core Reynolds number based on the free stream axial 

velocity may not be a very good estimate of the ratio of convection 

to diffustion transport. A more appropriate characteristic velocity 

would have been based on some maximum velocity, for example, 

~~2 2 

/a + V . Then instead of requiring 
Re I 'V. 0(1) 

for a numerical solution we would require 

Re h . 

IT ? 0(1) 

flO 

where Re is the core Reynolds number defined previously. This 
suggested that flows with high values of a and V could only 
be obtained with a smaller core Reynolds number. 

Several of the solutions obtained with Re = 100, and 
a = 1 are presented in Figures 19, 20, and 21. 

With V = 0.8944, the axial flow retardation produced 
the very slight stream surface expansion seen in Fig. 19a. This 
figure and Fig. 19b, which shows the variation of w on the axis, 
should be compared with Figs. 8a and 8g, which were obtained from 
the solution with identical a and V, but Re = 200. The 
larger shear forces, corresponding to the smaller Reynolds number, 
which develop when the axial flow near the axis is retarded, re- 
quire a stronger adverse pressure gradient to cause stagnation. 
Figures 20a and b show a well -developed breakdown obtained with 
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V = 1.095, Although higher initial swirl velocities are required, 
breakdown will still occur at this low Reynolds number. 

Figure 21a shows a breakdown obtained with subcritical 
upstream conditions in a very high swirl flow, with V = 1.342. 

The breakdown bubble is very large, and is followed by a low 
velocity region (Fig. 21b) in which a second reversal occurs. This 
rather large low velocity region may be an early stage in the 
development of a very large scale flow reversal similar to that 
obtained by Harvey (1962) in a vottex tube and mentioned in 
Section 1.2. 

4.6 Discussion 

Figure 22 shows the dependence of the minimum velocity 
on the axis, on V for all solutions obtained with Re = 200. 

Figure 11 displayed the case of a - 1 only. Solutions for which 
^ necessarily represent flows in which breakdown has occurred. 
The data points for equal values of a fall roughly into bands so 
that, in general, as a increases larger values of V are re- 
quired to yield a given minimum velocity. As pointed out in the 
previous section, this behavior is due to the dependence of the 
axial momentum of the flow near the axis, which must overcome an 
adverse pressure gradient, on a. The two a = 0.3 data points 
for V = 0.55 and 0.63 both represent solutions in which was 
a minimum at z = 0, and, therefore, equal to a. There was no 
retardation in these two cases since the large viscous farces at 
2 = 0 easily overcame the adverse pressure gradient developing 
downstream. 

Figure 23 is a plot of versus the maximum swirl 
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angle, <^, at z = 0 for the same solutions. As above, the data 
points fall into bands for equal values of a so that as a in- 
creases a given value .of corresponds to increasing values of 
(|). For convenience define ^ as the value of ^ at which = 0 
for each value of a. Then second order interpolation yields 
/ = 42.5°, 43.0°, 44.4°. and 47.2° for a = 0.3, 0.6, 1.0, and 

* 

1.4, respectively. The asymptotic nature of the variation of 
with a as a becomes small, which may be easily demonstrated by 
computing A4> /Aa for each of the three a intervals, suggests 

it * 

that there exists a minimum value of 4) such that for < (|> 
breakdown will not occur even for very small values of a. 

Figure 24 shows a point of versus 0-| . The collapse 
of the data, except for the two a = 0.3 cases which should be 
ignored since in those solutions no retardation at all took place, 
is particularly striking. For all values of a, the value of 
corresponding to = 0, which may be defined as 6-| , is very 
close to -0.05, e-j , therefore, gives a very good indication of 
the occurrence of breakdown. When 6-| < 0-j no breakdown occurs 
and, in contrast with the swirl angle, when e-j > B-j breakdown 
necessarily occurs regardless of the value of a. 

The appearance of a critical value of 0-j at which 
breakdown will occur is in agreement with the quasi-cylindrical 
theory of Mager (1972) discussed in Section 4.4, as is the relation- 
ships between and e-j implied by the data collapse. 

Figure 25 is a plot of versus 0-| for solution 
with Re = 100. The data again collapse onto a single curve. 

The amount of scatter in the data is less than in the previous 
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figure, perhaps since solutions at this lower Reynolds number will 
tend to be more accurate. It is clear from the figure that 
0^ « 0.01 for Re = 100, which suggests that e-j is a decreasing 
function of Reynolds number, in which case it could possibly 
approach Mager's 6^ for very high Reynolds numbers. This is, 
however, only a hypothesis which, due to the difficulty of obtain- 
ing solutions for high Reynolds number flows, cannot be substan- 
tiated. Also, as pointed out earlier, 0-j* is the value of e-| 
at Which large axial gradients will occur, and not necessarily the 
value for which reversed flow results. Therefore, even at high 
Reynolds numbers, e-| might be larger than e-j* and, thus, 
closer to the values predicted from these numerical solutions. 

Figures 26 and 27 show a comparison for several Reynolds 
numbers of w at r = 0 versus z for a = 1, and V = 0.63 
and 0.8944. The solutions at Re == 500 and 1000 were, as pre- 
viously mentioned, obtained using upwind differencing and are, 
at best, only suggestive of the actual flows. 

Figure 26 shows that the very mild retardation in the 
solution with Re = 200, a = 1, and V = 0,63, discussed in 
Section 4.2, becomes slightly more pronounced for decreasing Re 
and less pronounced as Re increases. This is in accord with 
the integral analyses of both Gartshore (1962,1963) and Mager 
(1972) In which Re enters in the solution only as a scale factor 
for z. Figure 27, however, shows that solutions which exhibit 
breakdown at small values of Re should continue to do so as Re 
Is increased. It is expected, therefore, that, while this study 
has obtained quantitative information for moderate Reynolds 
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numbers only, the conclusions reached here can be extended to 
higher Reynolds number flows. 

Regarding the applicability of Benjamin’s finite transi- 
tion theory, it must be noted that as mentioned before the solu- 
tions for Re = 200, a = 0.3, V = 0.8944; Re = 200, a = 0.6, 

V = 1.095; and Re = 100, a = 1, V = 1.342, which are dis- 
played in Figs. 13, 16, and 21, all exhibited breakdown with 
initial conditions which, based on Mager’s analysis, were sub- 
critical. The finite transition theory, however, requires that 
a vortex flow be supercritical before breakdown occurs, so that 
the numerical solutions suggest that the theory does not fully 
account for viscous vortex breakdown. 

In order to explain the appearance of an axisynmetric 
bubble followed by a spiral breakdown in highly swirling flows, 
and the appearance of only the spiral in flows with less swirl, 
Mager included the finite transition in his quasi-cylindrical 
analysis. He argued that the finite transition was a cross-over, 
ahead of the discontinuity, of the solution to a different branch 
of the solution curve with the same value of ©i ; and, that the 
discontinuity appears in the experiment as the spiral breakdown, 
while the cross-over, if it occurs, appears as the axisymmetric 
bubble. This argument, however, cannot explain the existence in 
Harvey's (1962) experiment of two axisymmetric breakdowns followed 
by a spiral in a steady flow, and the appearance, in Sarpkaya's 
(1972a) experiment, of similar phenomena in an unsteady flow. 

Based on the numerical solutions obtained in this study, 
and on the observations of Sarpkaya (1972a, b), the following 
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possible explanation of the various forms of vortex breakdown is 
proposed. 

In flows in which the swirl velocities are not large, 
the axial retardation, which results from diffusion of vorticity 
and occurs in all viscous vortex flows, may be such as to produce 
a region in the vortex core in which the flow is unstable to 
certain spiral disturbances. The introduction of such a disturb- 
ance would then result in the development of a spiral breakdown, 
Sarpkaya (1971b) has, In fact, suggested that the spiral is the 
result of an instability In flows with moderate swirl; and the 
agreement pointed out by Ludwieg (1965) between his stability 
theory (1962,1965) and the measurements of Kirkpatrick (1964) and 
Hummel (1965) upstream of the breakdown would then be explained. 
With high swirl, these instabilities might not occur and the 
retardation of the flow would then result in the axisymmetric 
bubble. The second retardation which occurs in the numerical 
solutions for large values of V may, however, result in a region 
of flow behind the bubble unstable to a spiral disturbance, in 
which case the flow could experimentally appear as an axisymmetric 
bubble followed by a spiral. Such flows have, of course, been 
observed, and since, based on the discussion of Section 4.2, 
there is no physical reason why more than two retardations cannot 
occur, the occurrence of a number of stable retardations, appear- 
ing as axisymmetric bubbles, followed by an unstable retardation, 
which would appear as a spiral, is completely in accord with this 
theory. 
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The spiral and axisyiranetric breakdowns may therefore 
only represent unstable and stable manifestations of the same 
physical phenomenon. 
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5.0 CONCLUSIONS 

Numerical solutions of the full Navier-Stokes equations have 
been obtained for vortex breakdown. The solutions show breakdown 
to be a necessary occurrence in unconfined vortex flows with a 
large enough value of a parameter 6^ which Is a measure of the 
flux of axial momentum deficiency, or flow force deficiency, in 
the vortex core. The value of 6-| In any flow is Independent of 
the subcrlticality or supercriticality of the assumed upstream 
conditions from which 6-| is determined, and solutions exhibiting 
breakdown were obtained with subcritical initial conditions. The 
finite transition theory which requires the flow upstream of 
breakdown to be supercritical, therefore, cannot fully account 
for vortex breakdown. 

The physical mechanism responsible for breakdown is the dif- 
fusion and convection of vorticity away from the vortex core, 
thereby requiring an increase in the pressure on the axis, and 
causing the retardation of axial flow. When the axial momentum 
of the flow near the axis is small compared to the pressure forces, 
which will be the case when e-j is large, the retardation is 
sufficient to produce stagnation and flow reversal and, thus, 
vortex breakdown. Along with pronounced axial retardation, the 
numerical solutions showed the breakdown bubble to be closed and 
contain very slowly rotating fluid. The maximum swirl velocity 
decreased very rapidly and the core radius increased suddenly in 
the breakdown region, which suggests that breakdown sould poten- 
tially be a useful means of lessening the danger aircraft wing-tip 
vortices present to following aircraft. 
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Due to the difficulty of obtaining numerical solutions at 
Reynolds numbers large enough for comparison, it was not possible 
to conclusively demonstrate that the failure of the quasi-cylindri- 
cal approximation, which occurs in the analysis of Mager (1972) 
for 0^ greater than a particular value 0^* independent of 
Reynolds number, always indicates the occurrence of breakdown. 
However, since the values of B-j at which breakdown occur tend 
to decrease toward as the Reynolds number of the solutions 

increases, and since the relationship between the minimum axial 
velocity and predicted by quasi-cyl indrical analysis has been 
verified in these numerical solutions, there is substantial evidence 
that the failure of the approximation, which is the result of large 
axial gradients, and which occurs in vortex flows exhibiting break- 
down, may also occur in those flows with substantial axial velocity 
retardation although slightly less than sufficient to produce 
stagnation. 

Finally, it has been suggested that breakdowns which appear 
in solutions with moderate values of swirl might be unstable to 
spiral disturbances, in which case the physical manifestation 
of the breakdown might be the single asymmetric spiral observed 
in the experiments of Harvey (1962) and Sarpkaya (1971a). Simi- 
larly, since solutions with large values of swirl exhibit a second 
axial flow retardation, if the first retardation were stable and 
the second unstable, the physical manifestation might be the axi- 
symmetric bubble followed by a spiral observed in the same experi- 


ments. 
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TABLE 1. REYNOLDS NUMBERS OF FINITE DOMAIN 
SOLUTIONS FOR a - V COMBINATIONS 
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TABLE 2. AND 0^ FOR SOLUTIONS WITH Re = 200, a = 1 


V 

“^max (degrees) 


0.63 

34.38 

-0.231 

0.80 

40.99 

-0.112 

0.85 

42.71 

-0.081 

0.8944 

44.17 

-0.056 

1.0 

47.36 

0. 

1.095 

49.95 

0.046 
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Fig. 8d. SWIRL VELOCITY PROFILES — Re = 
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AT FIXED RADIAL POSITIONS — R. 
6 ; (c) r = 0.714 ; (d 
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Fig. 16b. AXIAL VELOCITY ON AXIS vs. z — Re = 200, a = 0.6, V = 1.095 
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MINIMUM AXIAL VELOCITY AT AXIS VS,V 

Re = 200 


Fig. 22. 
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MINIMUM AXiAL VELOCITY AT AXIS 
VS. SWIRL ANGLE <#) 


g. 23. 
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Fig. 24. MINIMUM AXIAL VELOCITY AT AXIS VS. 0| 

Re = 200 
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APPENDIX A 

APPLICATION OF COORDINATE TRANSFORMATIONS 

To write equations (2.5) through (2.8) in terms of the vari- 
ables X and y, defined by the transformation equations (3.5) 
and (3.6), the following relations derived from the chain rule are 
required: 


-i. = 

9r 


3y dr 


(B-1) 


gr2 






(B-2) 


dx 

3x Hz' 


(B-3) 


2 2 2 

9 - 9 / dx .2 ^ 3 dS 

where, ft^om equations (3-5) and (3-6), 


(B-4) 


II 

e-ay 

ab 

4 = 

dr*^ 

e-2ay 

'1^ 

(B-5,6) 

dx 

H7 ■ 

-cx 

.2 

-2cx 


e 

cd 

II 

X Csi 
N 

o u 

e 

(B-7,8) 


Replacing derivatives in equations (2.5), (2.6), (2.7), and (2.8) 
by derivatives with respect to y and x using (B-1) through 
(B-4) , while letting 

f(y) = ^ , y(y) = ^ . s(x) = ^ , 


and t(x) = 

dz^ 


yields the transformed equations (3.7), (3.8), (3.9), and (3.10). 
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APPENDIX B 

FINITE DIFFERENCE FORMULATION 
DIFFERENCE EQUATIONS 

The ADI scheme applied to the momentum equations (3-7), (3-8), 
and (3-9), and the simple explicit scheme for the pressure equation 
(3-10) are explained as follows. The region 0£X<^1, 0£y£ 1/2 
in which a solution is to be obtained is overlaid with a grid system 
of N points in the x direction and M in the y direction. 

The location (x,y) of any grid point is given by the indices i ,j 
with x = (i-l)h, y = (j-l)k, where h and k are the uniform 
mesh width in the x and y directions (Fig. 3). u(x,y,t), 
v(x,y,t) and w(x,y,t) are represented by the functions 
and wj. defined at each grid point i,j. Superscript n is the 
number of time steps of length t which are equivalent to t. 
p(x,y,t) is represented by the function Pi+i/2,j+i/2 
a second grid staggered with respect to the first by a half mesh 
width in both coordinate directions. 

Centered differences are used for all spatial derivatives 
with the notation 





2h 
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where can be or w'Jj- These difference forms 

approximate the first and second derivatives in the x and the 
y directions, respectively. The pressure derivatives in the 
momentum equations at a point ij are approximated by an averaged 
differencing on the staggered grid so that 

Dqi Pij = 

* _ ^Pi+1/2..T-H/2-Pi>1/2.i-1/2^^^Pi-1/2..H1/2-Pj-1/2.j-1/2) 

°02 Pij = 2k 

represent 3p/3x and 9p/3y at i,j. Similarly, the velocity 
derivatives required to compute the pressure at a staggered grid 
point from the auxiliary continuity equation 1 + 1/2, j + 1/2 are 

approximated by the averaged differences 

* 

‘’m ”i+1/2.j+V2 

= 2h 

Do2(u/e)^^.^/2j+T/2 

= Tic 

With 6j = e(y), fj = f(y )9 9j = g(y)* = s(x), t^ = t(x), 

and using the definitions 


ij 

'ij 



u 


n+1/2 

ij 




the finite difference equation for the two ADI sweeps, with the 
pressure computation are 


First Sweep 




v?.e. 
1J 3 


1 r.2 


■ ^j‘’02Pij + l^fjD+2°-2*'ij ^ 






W tfK2°-2''ij ^ (9j+ejfj)Dg2V,.j 


+ s^+lD.,v|^ + t^Dg^vlj - v,j.e2] 


Wi 

+ fj“i jDo2«i j + 


1 r.2 


^i'^OlPij * Re f^o‘’+2°-2’'ij * j 


+ s^D+^D_^w!j + t^Dg^w^j] 


Pj+1/2.j+1/2 ~ Pi+1/2..i+1/2 

T/2 

l^®j+l/2^j+l/2°02^“^®^i+l/2,j+l/2 

* 

^i+l/2°0l”i+V2,j+1/2J^'* 
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Second Sweep 


u^.-u* 


+ f j^ijOo2‘'i j j - ''i j "j 

= - ^o‘'o2Pij " W [i^2°-2“ij " 

+ s?D+lD-l“i 3 + t^DoiUij - u^'j.e?] 


vijVv;. 

= ^ tf 3 -'>+ 2 °- 2 ''?j 




^ f93^®j^j^°02''ij 


wi'..-wj.. 


jDo 2 «^j + 


1 r ^2 




+ s^D^iD_^w!j + t^DQiW^j] 


Pi+1/2, 3+1/2 ■ Pi+l/2.j+1/2 
^F 72 

" ■ l^®j+l/2^j+1/2°02^“'''®^i+l/2,3+l/2 
®i+1/2°0l'*'i+1/2,3+l/2^^'^ 
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When rearranged, the finite difference momentum equations 
are of the form 


^^1-1 ^ ^•‘Ji ^ 


2 < i < N-1 


(B-1) 


Is obtained, therefore, by the inversion of a tridiagonal 
matrix. The criterion that such a matrix be well -conditioned, 

B- > - (A^+C. ), is always satisfied. Methods for inverting 
tridiagonal matrices when U-j and are specified are well 
known (Rlchtmyer and Morton, 1967). Extensions to other 
boundary conditions have been given by Roache (1972). 

The existence of two vectors and F. may be postulated 
such that 




(B-2) 


Combining (B-1) and (B-2) yields 


"i ■ B, + A. E.. 


''i 


1 


i ‘■i-1 ^i-1 




Comparison with (B-2) results in the recursion relations 


= 


B. + A,. E,.., 


-Ci 


(B-3),(B-4) 


A^. , B^, C^ and D^. are all known and application of the boundary 
condition at 1=1 yields and , as follows. 

Equation (B-2) at 1=1 is 
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U] = + E] ^2 (B-5) 

For a Dirichlet boundary condition, U-j = a, (B-5) requires 
F-| = a, E-| = 0, while for a Neumann boundary condition, 

3U/9y = 0, so to first order U 2 = and, therefore, (B-5) 
requires F-j = 0, E-| = 1 . 

With E-| and F^ specified, (B-3) and (B-4) determine all 
E^ and F^. , 1 < N. 

The boundary condition at 1 = N specifies as follows. 
For a Dirichlet condition = C where C is specified. 

For a Neumann condition, 3U/3y = 0, so U|^_-j = to first 
order and from 


“^N-l ° '"n-1 * ^N-1 ’^1 


(B-6) 


the requirement that 
F. 


U 


N 


N-1 

r^T 


N-l 


(B-7) 


is obtained. 

For a mixed condition, V + a(3V/3y) = 0, so that 


*^N + ) = 0 


(B-8) 


to first order, and therefore from (A-6) 




B - E 


N-l 


where 8 = 1^ 


With E^. , F^. and known, (B-2) may be used to compute 

, 1 < i < N-l . 
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Each sweep, therefore, proceeds as follows 

(1) calculation along each row or column of the 
coefficients A^, C. and D^. for each of 
the tridiagonal matrices for u*., and w.., 

(2) calculation of and and, hence, 
and F. for each variable from the boundary 
conditions at i = 1 (or j = 1), 

(3) application of the boundary conditions at i = N 
(or j = M) and calculation of u.., and 

for that row or column, 

(4) explicit pressure computation. 


